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Abstract 

In this paper, we develop a new approach to gravity-gradient noise subtraction for underground 
gravitational-wave detectors in homogeneous rock. The method is based on spatial harmonic expansions 
of seismic fields. It is shown that gravity-gradient noise produced by seismic fields from distant sources, 
stationary or non-stationary, can be calculated from seismic data measured locally at the test mass. Further- 
more, the formula is applied to seismic fields from stationary local sources. It is found that gravity gradients 
from these fields can be subtracted using local seismic measurements. The results are confirmed numeri- 
cally with a finite-element simulation. A new seismic-array design is proposed that provides the additional 
information about the seismic field required to ensure applicability of the approach to realistic scenarios 
even with inhomogeneous rock and non-stationary local sources. 
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Seismic waves produce perturbations of the gravity field, which are predicted to cause the 
so-called Newtonian or gravity-gradient noise (GGN) in gravitational-wave (GW) detectors 
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,\4\. Whereas the sensitivity of currently operating detectors is not limited by GGN 
151 16L |7|, Isl], second and third- generation detectors will be sensitive to gravity gradients at 10 Hz 
and below. Figure [U shows the current best estimate for some important noise contributions to 
the second-generation Advanced LIGO detector. The GGN curve is based on a model that ap- 
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FIG. 1: Current estimates for some important noise contributions to Advanced LIGO (generated with the 
GWINC matlab script). 

proximates a typical spectrum of seismic surface waves measured at the Hanford site of the LIGO 
detectors, and is based on the characteristics of local geology. Third-generation GW detectors 
will be designed with enhanced sensitivity below 10 Hz based on improved sus pen sion systems 
and optimized material properties to mitigate the seismic and thermal noise Ba, uM. Moreover, 
quantum-non-demolition techniques are being investigated to cancel part of the optical quantum 
noise llll [l2[ [isll . This leaves the question whether the GGN, which is directly imprinted on the 
test-mass motion without the possibility to build an isolation system, can be decreased. One ob- 
vious improvement would be to identify a detector site with a comparatively low level of seismic 
noise, which also includes the possibility to construct the detector under ground. Underground 
seismic noise at depths of about 1 km is known to be an order of magnitude weaker than surface 
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but further GGN mitigation by two orders of magnitude is required 
to achieve good sensitivities at frequencies close to 1 Hz. The solution is to subtract from the GW 
data an estimate of GGN based on seismic measurements. So far, the assumption was that a 3D 



array consisting of several hundred seismometers deployed around each test mass extending over 



several seismic wavelengths would be required lll7h . In this paper, we will show for homogeneous 



media that GGN produced by seismic fields from distant sources and from stationary spherical 
waves can be subtracted using relatively few seismometers positioned at the test mass. 

The calculation is based on a plane- wave expansion of the seismic displacement field ^(r, t). 
Let us first consider the simple case of a plane pressure wave (P-wave) from a distant source, in 
which case its frequency u and wave number k = \k\ are related by lu = kcp, Cp being the speed 
of pressure waves. The P-wave with initial phase 0o at the origin r = produces longitudinal 
displacement along the direction Ck = k/k: 

^(f, t) = Ckil exp(i((/.o + ujt-k-r)) (1) 

Assuming that the test mass is located at the origin, the lowest order perturbation of the gravity 



field can be calculated using the dipole formula [Il7h 



5a{t) = Gpo j (e(f, t) - 3(6-; ■ f(r, , (2) 

where G is Newton's constant, and po denotes the mean density of the rock. The integral is carried 
out over the entire volume of the medium. Inserting the displacement field of the P-wave, it is 
straight-forward to solve the integral in spherical coordinates. If one allows for an arbitrary range 
of integration in radial direction, then the solution is found to be 
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The inner radius ri can be interpreted as the radius of a spherical cavity surface. A finite outer 
radius r2 is needed to evaluate corrections of the result when the integration volume is bounded. 
The integral as a function of distance is plotted in Fig. |2l Instead of considering gravity gradients 
as a function of the radius r2 of a spherical rock volume, one could ask for the convergence curve 
of gravity gradients in an infinite rock volume. The difference is that contributions of gravity- 
gradients from the outer surface at r = r2 need to be subtracted from Eq. dS]). For plane pressure 
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FIG. 2: The figure displays a plot of the normalized gravity-gradient integral s{kr2) — s(0) which appears 
in Eq. and the convergence curve of gravity gradients from plane waves in an infinite rock volume 
as function of integration distance. The two functions differ by the surface term Eq. dH) which needs to be 
subtracted from Eq. Q to obtain the convergence curve. This example demonstrates that in certain problems 
one needs to take care of surface terms when using the dipole result Eq. Q. Numerical results obtained 
from a finite-element simulation agree well with the analytic result Eq. Here, the simulated seismic 
field was constructed by adding 50 plane P-waves with random directions of propagation. More details can 



be found in um- 

waves, the surface contribution reads 

5CfW =47rGpor(0,t) 
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One does not have to subtract the surface contribution from Eq. ([3]) when comparing analytical 
and numerical results, but it becomes necessary in reality when evaluating gravity-gradient contri- 
butions as a function of distance. 

The gravity-gradient integral for a plane P-wave, with or without surface term, is of the form 
^^(0, t) ■ f{k). Therefore, it is possible to add multiple plane P-waves with arbitrary directions 
of propagation to construct any stationary P-wave field for a specific wave number k, and the 
vector ^^(0, t) becomes the sum of displacements of all these waves at the origin. So far, we 
have not made use of the fact that the field Eq. ([T]) obeys the plane-wave equation. The gravity- 
gradient formula Eq. ([3]) can be read as an identity in fc-space, interpreting k as an independent 
spatial-harmonics variable, and ^(0, t), 6a{t) as amplitudes ^( A;, t), a{k, t) of the respective spatial 
spectrum. This makes it valid for arbitrary stationary or non- stationary fields from local or distant 



sources. We will come back to this in the next paragraph. However, the usefulness of the formula 
for the problem of GGN subtraction is most obvious for fields from distant sources that obey 
the plane- wave equation enforcing uj = kcp. Translating wave numbers into frequencies, it is 
clear that the function f{k = uj/cp) can be multiplied in frequency domain after calculating the 
temporal Fourier transform of Eq. ([3]). In other words, it is possible with local data and simple 
transformations in frequency domain to calculate GGN for a certain range of frequencies. In 
practice, the applicability is limited to frequencies / > Cp/H, H being the depth of the test mass, 
since surface effects become important at lower frequencies. A plane-surface correction to Eq. ([3]) 
is required to improve low-frequency predictions of GGN. 

A concise form of Eq. ([3]) can still be maintained if shear waves, which produce GGN ex- 
clusively through surface displacement, are added to the total displacement ^{r, t) = ^^{r,t) + 
^^{r, t). For an arbitrary spatial spectrum, one obtains 



In practice, disentangling the shear and pressure modes is accomplished by a small array of sensors 
around the test mass, or by instruments which specifically respond to pressure or shear waves 
(e.g. strainmeters and rotational instruments, or a measurement of the test-mass position relative 
to the cavity walls). Once the modes are identified, the integral does not explicitly depend on the 
direction Ck of propagation. In this way, integration over different directions can be carried out, 
leaving the integral over the wave number k 
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in terms of the "direction-averaged" amplitudes: 
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where the volume element is written d^k = k'^dkdfl. As for the case of a single plane P-wave, 
Eq. Q is easy to evaluate in frequency domain for fields satisfying the plane- wave equation. This 
result is important, not only because it greatly simplifies the task of subtracting GGN produced 
by seismic fields from distant sources in homogeneous rock, but also because it serves as starting 
point to investigate the subtraction problem for fields from local sources in terms of their spatial 
spectrum. For the discussion of seismic fields from local sources, it will be useful to derive an 
approximation to the gravity-gradient formula: 
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which neglects terms of order 0{{kri)^), 0{{kr2Y) in the square brackets. The idea is to approx- 
imate the cavity as small kri <^ 1 and the outer surface as large kr2 ^ 1. In the limit ri ^ and 
r2 — oo, the formula reduces to 

4-7r -» 

<5a(t) = — Gpo(2r(0,t)-e'(0,t)) (9) 

The latter equation holds for any field in homogeneous media under the assumption of infinitely 
small cavity surface and infinitely large outer surface. In other words, whenever the lowest order 
corrections in Eq. ([8]) are negligible, then its limit Eq. (HI) can be applied which is a simple function 
of local seismic data. Whether the approximation is valid depends on the wave number and the 
corresponding direction- averaged amplitude. If the spatial spectrum of a field contains non-zero 
amplitudes ^d.a.(^, t) for all k, then the question is whether the amplitudes are negligible at high 
wave numbers kr^ > 1 and small wave numbers fcro < 1. 

As an extreme example for time- variable distant sources, one can test the applicability of Eq. ^ 
with seismic fields from earthquakes. Comparing numerical integrations of Eq. (O with Eq. Q in 
the case of a real earthquake signal, we found a better than factor 20 reduction throughout the 
entire waveform. For the East- West displacement and a small stretch of the waveform, the result 
can be seen in Fig.[3l As one would expect, the error is largest at times when the earthquake signal 
changes fastest. The subtraction is improved to the level of numerical accuracy, if a frequency- 
dependent subtraction is applied based on Eq. 

Our next step is to show that Eq. ^ can also be applied to calculate gravity-gradients gener- 
ated by seismic fields from stationary local sources. Here one has to take into account that the 
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FIG. 3: The plot shows part of the simulated GGN from an earthquake measured on February 18th, 2009 
at the Kamioka station GIFHIO in Japan. The earthquake had a magnitude of M = 5.2. Its epicenter was 
at a distance of 124 km to GIFHIO. Subtracting the GGN estimate Eq. Q from a numerical integration of 
Eq. (O, a factor 20 reduction of GGN is achieved. 

plane- wave relation u = kc does not hold for spherical waves. The wave number k needs to be 
considered as an independent parameter used in the spatial harmonic expansion of the field. As an 
example, the stationary P-wave fields from local sources is investigated using the spherical wave 
originating from f=fo: 
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with uj = koCp, ko{r — tq) = ko\r — ro\, and Cr = (r — ro)/|r — ro|. Typical local underground 
sources of seismic waves include Rayleigh-scattering centers, pumping stations, and ventilation 
systems that produce fields with radiation patterns of varying complexity. Therefore, other types 
of fields from local sources need to be studied in the future, but here our main intention is to 
demonstrate how the spatial harmonic expansion can be applied to more general fields. For the 
chosen spherical wave, one expects a spatial spectrum that extends over a range of wave numbers 
k, with spatial amplitudes peaked at /c = k^. A straight-forward calculation leads to the following 
spatial Fourier transform of Eq. (flOl) . valid for k ^ k^: 



e{k,t) = 2TT^oe 
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The corresponding direction- averaged amplitudes are obtained by inserting the spatial amplitudes 



7 



350 




Wave Number [k/k 



FIG. 4: The figure shows a plot of the modulus |^d.a. I of the direction-averaged amplitudes of a spherical 
P-wave. The distance to the center of the spherical wave is tq = 1, and the displacement amplitude is 
normalized to = 1- 
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The spectrum is plotted in Fig. HI The plot shows that the spectrum is indeed peaked atk = ko and 
that the amplitudes fall rapidly towards higher and lower wave numbers: 
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In the following, we test whether Eq. ^ produces accurate predictions of GGN from spherical 
waves. As an example, let us assume that the cavity has a radius ri = 20 m and that the test mass is 
2 km underground, so that Eq. ([8]) can be evaluated choosing r2 = 2 km. Now, let the wave number 
of the spherical wave he ko = 27r/(l km), which corresponds to a 3 Hz wave propagating in rock 
with P-wave speed 3km/s. We find that (fcori)^ ^ 0.016 and {kor2Y ~ 160. Therefore, at A; = ko, 
both corrections in Eq. ([8]) can be neglected given a subtraction goal of two orders of magnitude. 
Whether amplitudes at higher and smaller wave numbers are negligible, in the sense that Eq. (HI) 
is still a valid approximation of the exact solution Eq. is determined by evaluating the integral 
Eq. ^ over wave numbers k where the approximation Eq. ([8]) does not hold. Alternatively, using 
a finite-element simulation, one can directly compare gravity gradients from Eq. Q and Eq. 



8 



One finds that the relative deviation depends on the distance of the centers of the spherical waves 
to the outer surface. The closer the center is to the surface, the greater is the deviation, since \ow-k 
surface corrections that are explicitly neglected in Eq. ^ become more significant relative to the 
weak gravity gradient produced by more distant sources. The simulation showed that the relative 
deviation is smaller than 0.02 for all spherical waves that originate at least one wavelength away 
from the outer surface, and increases to 0. 1 for spherical waves which originate very close to the 
outer surface. This observation may be relevant in practice for strong surface sources. For this 
reason, a surface array of seismometers is required to guarantee sufficient subtraction of gravity 
gradients from individual strong surface sources. 

In this paper, an approximate time-domain estimate of GGN, Eq. and a frequency-domain 
estimate based on Eq. Q, that is exact for seismic fields from distant sources, have been developed 
for homogeneous media. If ideal conditions as homogeneity of the rock and stationarity of local 
sources do not apply, then there are several options to adapt the seismic array and to gather the 
information needed to improve the gravity-gradient model for realistic seismic fields. In addition 
to the seismometers at the test mass, sensors can be deployed near large geo logic faults and other 
inhomogeneities that produce significant scattering of the seismic field [ISQ. A surface array of 
seismometers would be necessary to subtract gravity gradients at very low frequencies where the 
Earth surface cannot be treated as far away from the test mass, and to improve subtraction of strong 
local surface sources. A dense seismic array around the test mass could be used to measure the 
spatial amplitudes of the seismic field as a function of time, which would allow us to evaluate 
Eq. Q in a more accurate form jiol . and improve subtraction of all local sources. 

We gratefully acknowledge the support of LIGO, and thank the operators of the Kiban-Kyoshin 
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